%%Test Program for Strang RC Circuit Using ODE45
% Implements Strangs modified nodal analysis with:
% Extension to general reactance (i.e. C and L) and:
% Uses matlab ODE solver
% NOTE: Term KVL is not quite right in Strang
% He doesn't really use the voltage drop around loops,
% Only KCL. KVL is used when the edge-loop matrix (i.e. the curl)
% is used, as in more traditional, ECE type circuit analysis (!)
% Also, note that the matrixes G,C and L are singular (i.e. they have
% rows with zero entries. This is not a problem since they only occur in the form
% A'GA, etc.-- they are projectors for the relevant edge variables (!)

%% Reference Material
% http://cns.bu.edu/~eric/readings/strang.pdf

%% Circuit Diagram
% http://cns.bu.edu/~eric/readings/RC.jpg

%%  Variables
% The following matrixes and vectors are needed
%
% Sparse matrixes are used throughout the code
%
% * A Incidence matrix  node by edge
% * A_0 Reduced incidence matrix nodes by edge-1
% * b Voltage source vector b edge by 1
% * f Current source vector node by 1
% * G Conductance matrix G edge by edge
% * C Capacitance matrix edge by edge
% * y Current vector  node by 1
% * x Potentials node by 1
% * e Potential differences edge by 1

%% Kirchoff KCL and KVL
% KCL and KVL may be succinctly stated in matrix vector form as  follows.
%
% Ohms law (constitutive) appears below in conductance form: 

%% Ohms Law
% $$y_r=Ge$$
%
%% A capacitive constituitive law
%%
%
% $$y_c=\mbox{-}C\dot{e}~=~\mbox{-}C\frac{d}{dt}(b~\mbox{-}A_0 x)~=~ CA_0\dot{x}$$
%

%% A capacitance matrix C allows reactive circuits to be handled.
%
% The ohmic current, together with the displacement (capacitive) current, simply
% balances the current sources:
%%
%%
%
% $$A_0^tG(b~ \mbox{-}~A_0x)+A_0^t CA_0\dot{x} ~=~ f$$
%
%% SUMMARY OF KCL/KVL/OHM

%% Edge voltage drops
% $$e=b~ \mbox{-}~ A_0x$$

%% Capacitive (displacement current) 
%
% $$y_c~=\mbox{-}~C\dot{e}=CA_0\dot{x}$$ 

%% Ohm's Law
% $$y_r=Ge$$

%% KCL
%
% $$A_0^ty=f$$

%% KCL/KVL with only resistive edge variables according to Strang
%%
% $$A_0^tG(b~\mbox{-}~A_0x)=f$$
%

%% KCL/KVL with resistive and capacitive  elements according to ELS
%%
%
% $$A_0^tG(b~\mbox{-}~A_0x)-A_0^t CA_0\dot{x} ~=~ f$$ %- sign on
% displacement current

% ELS 11/1/2006

% ELS 9/14/2007 added in ode45 code and
% could add inductance by appropriate constitutive law:
%% Outline extension to inductive circuits
%
% $$\dot{y}=Le$$
%
% This makes the system second order, but it is already in
% the "dimension reduced form of y,\dot(y) for first order solution!
% The full equation including inductors can be written as:
%
% $$ A^t G (b - A_0x) + A^tCA_0 \dot{x} + A^t \int^t L A_0 x ~=~ f$$
%
%
% $$ y_r + y_C + y_L $$
%
% To get this into second order ODE form, take time derivative
%
% $$ A^t G (b - A_0 \dot{x}) + A^t C A_0 \dot{\dot{x}} + A^t L A_0 x ~=~ f $$
%
% Let z = stacked of x, and dx/dt, then
% 
% $$\dot{z}= M Z + B$$
%
%%Incidence matrix
% A= -1 1;1 -1
% A_0= [-1 1]^t
i=[];j=[];s=[]; %initialize to empty
%
i=[i 1]; j=[j 1]; s=[s -1];
i=[i 1]; j=[j 2]; s=[s 1];
i=[i 2]; j=[j 1]; s=[s 1];
i=[i 2]; j=[j 2]; s=[s -1];
A=sparse(i,j,s);
%
%% Reduced incidence matrix: chop last column
A_0=A(:,1); % chop last column to form reduced Incidence matrix
%
%% Conductivity matrix and Capacitance Matrix
R1=10; %10 Ohms
R2=5; %5 Ohms
C=.01; %.01 Farad
G=[1/R1 0;0 0];
C=[0 0; 0 C];
G=sparse(G);
C=sparse(C);
%% Voltage sources
%Battery V_0 pos up
%two Resistors R_1=10 and R_2=5
%
V_0=10; % 10 volt source
b=[-V_0;0]; % b is a column vector
b=sparse(b);
%% Solve using brute force Euler
%
% $$A_0^t G (b~\mbox{-}~A_0x) - A^tCA_0\dot{x}~-~f ~=~ 0$$
% 
x=zeros(100,1);
t=zeros(100,1);
x(1)=0; t(1)=0;
delt=0.01;
for j=2:length(x)
x(j)=x(j-1)+delt*A_0'*G*(b-A_0*x(j-1))/(A_0'*C*A_0); %note fantastically busy loop!
t(j)=t(j-1)+delt;
end
figure,plot(t,x); title('Euler by Hand'); axis([0 1 0 10]);ylabel('Volts');xlabel('Time-Secs');

%% Set up for ODE 45
%
% $$A_0^t G (b~\mbox{-}~A_0x) \mbox{-} A^tCA_0\dot{x}~\mbox{-}~f ~=~ 0$$
%
% $$\dot{x}={(A_0^t C A_0)}^{\mbox{-1}}$$
%
% $$A_0^tG(b~\mbox{-}~A_0x)~\mbox{-}~{(A_0^t C A_0)}^{\mbox{-1}}f$$
%
% $$INV={(A_0^t C A_0)}^{\mbox{-1}}$$
%
% $$CC=INV*A_0^tGb - INV*f$$
%
% $$DD=INV*A_0'*G*A_0$$
%
% $$\dot{x}= CC \mbox{-} DD*x$$
%

%% Set up vectors x,b,f,e and summary matrixes CC,DD, and INV
% Rank of A_0 is node x edge
% so, Rank of node variables (y,x) is second dimension of A_0
% Rank of edge varialbes (b,e) is first dimension of A_0
%
[r,c]=size(A_0);
%If A_0 has c cols, then node vars need c rows
y=zeros(c,1); y=sparse(y) ; %current node vector
f=zeros(c,1); f=sparse(f); %current source node vector
x=zeros(c,1); x=sparse(x); %potential node vector
x_0=zeros(c,1); x_0=sparse(x_0); %initial values node vector
%
INV=inv(A_0'*C*A_0);
CC=INV*A_0'*G*b - INV*f;
DD=INV*A_0'*G*A_0;
gfunc = @(t,x) CC - DD*x;
[t,x]=ode45(gfunc,[0:0.01:1],x_0(1:c));
figure,plot(t,x); title('ODE Solver'); axis([0 1 0 10]);ylabel('Volts');xlabel('Time-Secs');


